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1 .  INTRODUCTION 

i 

Me  consider  singularly-perturbed  two-point  boundary  value  problems  for 
nonlinear  vector  systems  of  the  form 

x  -  f(x,y,t,e)  ,  ey  -  g(x,y,t,e)  ,  0  <  t  <  1  (la,b) 

a(x(0)*,y(0),e)  -  0  ,  b(x( 1) ,y( 1) ,e)  =  0  (lc,d) 

where  x,  y,  a,  and  b  are  vectors  of  dimension  m,  n,  q,  and  r  =  m  +  n  -  q, 
respectively.  We  seek  to  find  limiting  solutions  of  problem  (1)  as  the  small 
positive  parameter  e  tends  to  zero;  however,  to  do  this  in  complete  generality 
is  very  difficult  and  beyond  the  grasp  of  our  current  understanding.  Thus,  we 
simplify  problem  (1)  considerably  by  assuming,  in  addition  to  natural  smooth¬ 
ness  hypotheses,  that  (1)  g,  a,  and  b  are  linear  functions  of  the  fast 
variable  y,  i.e. 

g(x.y.t.e)  -  gi(x,t,e)  +  G2(x,t,e)y  (2a) 

a(x(0) ,y(0) ,e)  »ai(x(0),e)  +  A2(x(0) ,e)y(0)  (2b) 

b(x( 1) ,y(l) ,e)  -  b1(x(l),e)  +  B2(x( 1) ,e)y(l)  (2c) 

(il)  that  G2(x,t,e)  has  a  hyperbolic  splitting  with  k  >  0  stable  eigenvalues 
and  n  -  k  >  0  unstable  eigenvalues  for  all  x  and  0  <  t  <  1,  and  (ill)  that  q  > 
k  and  r  >  n  -  k. 

With  the  assumed  hyperbolic  splitting,  we  wjuld  expect  y  to  vary  rapidly 
relative  to  (the  slow  vector)  x  in  narrow  boundary  layer  regions  near  both  t  » 
0  and  1.  We  thus  seek  limiting  solutions  having  the  form 

x(t,e)  -  X(t)  +  0(e)  ,  y(t,e)  -  Y(t)  +  u(t)  +  v(a)  +  0(e)  (3a, b) 

where  the  initial  layer  correction  y(t)  and  the  terminal  layer  correction 
v(o),  respectively,  decay  to  zero  as  the  stretched  variable 

t  -  t/e  or  a  -  (l-t)/e  (4a, b) 
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tend  to  infinity.  The  limiting . solution  X(t),  Y(t)  within  0  <  t  <  1  must 
necessarily  satisfy  the  reduced  system 

X  =  f (X,Y,t,0)  ,  0  =  g(X,Y,t ,0)  (5a, h) 

Because  G2  is  everywhere  nonsingular,  we  can  use  Eqs.  (2a)  and  (5b)  to 
determine 

Y(t)  =  -G2-1U,t,0)gl(X,t,0)  (6) 

in  a  locally  unique  way,  and  there  remains  the  mth  order  nonlinear  differen¬ 
tial  system  (Eq.  (5a))  for  determining  X(t). 

In  order  to  completely  specify  the  reduced  solution  we  must  prescribe  m 
boundary  conditions  for  equation  (5a).  We  do  this  by  providing  a 
''cancellation  law"  which  selects  a  combination  of  q-k  initial  conditions  (Eq. 
(2b))  and  of  r  -  n  +  k  terminal  conditions  (Eq.  (2c))  to  be  satisfied  by  X  and 
Y.  In  Section  2  we  present  a  numerical  procedure  for  determining  the  boundary 
conditions  for  the  reduced  system  that  uses  an  orthogonal  matrix  E(x,t)  to 
reduce  the  matrix  G2(X(t),t,0)  to  a  block  tridiagonal  form  so  that  the  stable 
and  unstable  elgenspaces  may  be  separated.  The  boundary  layer  corrections 
M(t)  and  v(o)  in  Eqs.  (3)  compensate  for  the  cancelled  initial  and  terminal 
conditions,  respectively,  and  they  can  be  determined  once  X(t)  has  been 
computed  (cf.  Section  2).  This  process  avoids  complicated  matching 
procedures. 

In  Section  3  we  discuss  a  numerical  procedure  for  determining  the 
asymptotic  approximation  (Eq.  (3))  which  uses  the  general  purpose  two-point 
boundary  value  code  COLSYS  to  solve  the  reduced  problem  and  then  adds  numeri¬ 
cal  approximations  to  the  boundary  layer  corrections.  This  approximation  is 
considerably  less  expensive  to  obtain  than  solving  the  full  stiff  problem 
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numerically  and  It  has  the  advantage  of  Improving  In  accuracy,  without  any 
additional  computational  cost,  as  the  small  parameter  e  tends  to  zero.  How¬ 
ever,  when  e  is  only  moderately  small  our  asymptotic  approximation  may  not  be 
sufficiently  accurate  for  some  purposes,  so  we  have  developed  a  procedure  (cf. 
Section  3)  that  generates  an  improved  solution  by  using  COLSYS  to  solve  the 
complete  problem  (Eqs.  (1)  and  (2))  with  our  asymptotic  approximation  as  an 
initial  guess.  In  order  for  this  approach  to  succeed  we  must  also  provide 
COLSYS  with  an  initial  nonuniform  mesh  that  is  appropriately  graded  in  the 
boundary  layers  (cf.  Ascher  and  Weiss*)  and  we  give  an  algorithm  for 
constructing  such  a  mesh  in  Section  3.  While  our  procedure  does  not  appear  to 
be  optimal,  we  show  by  an  example  Involving  the  deformation  of  a  nonlinear 
elastic  beam  (cf.  Section  4)  that  it  does  offer  some  advantage  over  the  more 
standard  approach  of  continuation  in  e,  where  one  starts  with  a  large  value  of 
e  (e.g.  £  -  1)  and  a  crude  initial  guess  and  reduces  e  in  steps  so  that  the 
mesh  is  gradually  concentrated  into  boundary  layer  regions. 

We  close  Section  4  with  a  second  nonlinear  beam  example  that  is  beyond 
the  capabilities  of  our  present  methods  because  the  matrix  G2  is  a  function  of 
y.  Flaherty  and  O'Malley^  analyzed  this  problem  and  showed  that  its  solution 
becomes  unbounded  as  e  ♦  0.  We  include  the  numerical  solution  of  this  problem 
in  this  paper  in  order  to  show  one  of  the  many  challenging  effects  that  can 

*U.  Ascher  and  R.  Weiss,  "Collocation  For  Singular  Perturbation  Problems  I: 
First  Order  Systems  With  Constant  Coefficients,”  Technical  Report  81-2,  Dept. 
Comp.  Scl.,  University  of  British  Columbia,  1981. 

2j.  E.  Flaherty  and  R.  E.  O'Malley,  Jr.,  "Singularly-Perturbed  Boundary  Value 
Problems  For  Nonlinear  Systems,  Including  a  Challenging  Problem  For  a  Non¬ 
linear  Beam,"  Proceedings,  Conference  on  Slngulare  Storungstheorle  mlt 
Anwendungen,  Oberwolfach,  1981.  ~ 


occur  with  singularly-perturbed,  problems. 

Finally,  in  Section  5  we  discuss  our  results  and  present  some  suggestions 
for  future  investigations. 

2.  ASYMPTOTIC  APPROXIMATION 

In  order  to  calculate  the  boundary  conditions  for  the  reduced  problem 
(Eqs.  (5a)  and  (6))  and  the  boundary  layer  corrections  y(i)  and  v(o)  we 
calculate  the  Schur  decomposition  of  the  matrix  G2  at  t  »  0  and  t  =  1.  In 
particular,  at  t  =  0  we  find  an  orthogonal  matrix  E(x(0))  such  that 

T_(x(0) )  U(x(0) ) 

G2(x(0) ,0,0)E(x(0) )  =  E(x(0) )  (7) 

0  T+(x(0))  • 

where  T_  is  k  x  k  and  upper  triangular  with  the  stable  eigenvalues  of  G2,  and 
T+  is  upper  triangular  with  the  n-k  unstable  eigenvalues  of  G2.  The 
decomposition  (Eq.  (7))  can  often  be  obtained  analytically;  however,  when  this 
is  not  possible  or  practical  it  can  be  determined  numerically  by  using  the  QR 
algorithm  (cf.  Golub  and  Wilkinson^  and  Ruhe^  for  specific  procedures). 

We  partition  E  after  its  kth  column  as 

E  *  [ E_  E_]  (!) 

and  note  that  E_  spans  the  stable  eigenspace  of  G2  at  t  =  0  and 

P  =  E_  E-T  (9) 

^G.  H.  Golub  and  J.  H.  Wilkinson,  "Ill-Conditioned  Eigensystems  and  the 
Computation  of  the  Jordon  Canonical  Form,"  SIAM  Review  18  (1976),  pp. 

578-619. 

^A.  Ruhe,  "An  Algorithm  for  Numerical  Determination  of  the  Structure  of  a 
General  Matrix,"  BIT,  10  (1970),  pp.  196-216. 


is  a  projection  onto  this  eigenspace. 

Near  t  »  0,  we  assume  that  the  terminal  layer  correction  v  is  negligible, 
substitute  the  asymptotic  approximation  (Eq.  (3))  into  the  differential 
equations  (Eqs.  (la,b)),  use  the  reduced  system  (Eq.  (5)),  and  retain  only  the 
leading  order  terms  to  find  that  y(r)  satisfies  the  conditionally  stable 
system 

du/di  -  G2(0) u  (10) 

where  (here  and  below)  we  use  the  argument  t  to  denote  conditions  evaluated  at 
x(t)  -  X(t),  t,  and  e  »  0,  e.g., 

G2(0)  :»  G2(X(0),0,0)  (11) 

Integrating  Eq.  (10) 

G2(0) t 

u(t)  -  e  w(0)  (12) 

We  require  that  u(r)  decays  as  t  increases  and  this  will  be  the  case  provided 
that  u(0)  is  in  the  stable  eigenspace  of  G2(0) ;  thus,  using  Eq.  (9)  we  require 

U(0)  -  P(0)w(0)  -  L(0)Lt(0)v(0)  (13) 

Using  Eqs.  (3),  (13),  and  (2b)  in  Eq.  (lb)  we  find  that  the  limiting 
initial  conditions  have  the  form 

«l(0)  +  A2(0)  [Y(0)  +  E-(0)E-T(0)jj(0)J  -  0  (14) 

We  assume  that  A2(0)E_(0)  has  its  maximal  rank  k  and  construct  a  q  x  q  matrix 

L*  -  {L_t  L_t)  (15a) 

that  reduces  it  to  row  echelon  form,  i.e., 

V_ 

A2(0)E-(0)  -  (15b) 

0 
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where  V-  is  k  x  k  and  nonsingular.  Multiplying  Eq.  (14)  by  L  and  using  Eqs. 
(13)  and  (15)  gives  the  initial  layer  jump  and  the  q-k  initial  conditions  for 
the  reduced  problem,  respectively,  as 

p(0)  =  -E_(0)V__1L_[ai(X(0) ,0)  +  A2(X(0) ,0)Y(0) ]  (16a) 

and 

$(X(0))  L_[ai(X(0),0)  +  A2(X(0) ,0)Y(0) ]  =  0  .  (16b) 

We  find  the  terminal  layer  jump  and  the  r  -  (n-k)  terminal  conditions  for 
the  reduced  problem  in  an  analogous  fashion  with  the  exception  that  we  define 
E(x(l))  such  that  __ 

A  A 

T+(x( 1) )  U(x( 1)) 

G2(x( 1) , l,0)E(x( 1) )  =  E(x(l))  -  (17) 

0  T_(x( 1) ) 

which  we  partition  after  its  (n-k)th  column  as 

E  -  [Ef  M  (18) 

A  A 

In  parallel  with  Eqs.  (7)  and  (3),  the  matrices  T_,  T+,  and  E+  contain  the  k 
stable  eigenvalues,  the  n-k  unstable  eigenvalues,  and  span  the  unstable  eigen- 
space,  respectively,  of  G2  at  t  =  1.  Our  reasons  for  switching  the  positions 
of  the  matrices  containing  the  stable  and  unstable  eigenvalues  of  G2  is  that 
there  is  no  simple  and  stable  computational  procedure  for  finding  a  set  of 
vectors  that  span  a  given  subspace  and  are  not  in  the  leading  columns  of  in 
orthogonal  matrix  like  E  (cf.  Golub  and  Wilkinson-*). 

^G.  H.  Golub  and  J.  H.  Wilkinson,  "Ill-Conditioned  Eigensysteras  and  the 
Computation  of  the  Jordon  Canonical  Form,”  SIAM  Review  18  (1976),  pp. 

578-619. 
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Now,  following  the  procedure  that  we  used  for  the  initial  layer,  we  find 


that  the  terminal  layer  correction  satisfies 

G2(1)o 

v(o)  *  e  v(0)  ( 19) 

In  order  for  v(o)  to  decay  as  a  increases  we  require  v(o)  to  be  in  the 
unstable  eigenspace  of  G2(l);  thus,  we  take 

v(0)  -  Q(l)v<0)  -  E+(l)K+-TOM0)  (20) 

where  Q  is  a  projection  onto  the  (n-k)  dimensional  unstable  eigenspace  of 
G2( 1) - 

We  assume  that  B2(1)E+(1)  has  its  maximal  rank  n-k  and  find  a  r  x  r 
matrix 

rT  -  [R+T  r+T]  (21a) 

that  reduces  it  to  the  row  echelon  form 


R+ 

B2( l)E+( 1)  - 

V+ 

Rf 

0 

where  V+  is  (n-k)  x  (n-k)  and  nonsingular.  Multiplying  Eq.  (Id)  by  R,  using 
Eqs.  (2c),  (3),  (20),  and  (21),  and  retaining  only  the  leading  order  terras  we 
find  the  terminal  layer  jump  and  the  r  -(n-k)  terminal  conditions  for  the 
reduced  problem,  respectively,  as 

v(0)  -  -E+(l)V+_1Rf[bl(X(l),0)  +  B2(X( 1 ) ,0)Y( 1 ) ]  (22a) 

and 

*(X(1))  R+[bi(X(l),0)  +  B2(X( 1) ,0)Y( 1) ]  -  0  (22b) 


in  Che  interest  of  brevity,  we  have  omitted  several  details  of  our 
construction  and  have  not  attempted  to  justify  the  asymptotic  validity  of  our 
procedure.  These  topics  will  be  the  subject  of  a  forthcoming  paper  by 
O'Malley  and  Flaherty. ^ 

3.  NUMERICAL  PROCEDURE 

Our  computational  procedure  consists  of  first  solving  the  reduced  problem 
(cf.  Eqs.  (5a),  (6),  (  16b)  ,  and  (22b))  numerically  and  then  adding  any 
boundary  layer  corrections.  Since  the  reduced  problem  is  not  stiff  we  can  use 
any  good  code  for  two-point  boundary  value  problems  (cf.  Childs  et  al.^)  and 
we  have  chosen  '"o  use  the  collocation  code  COLSYS  of  Ascher,  Christiansen,  and 
Russell . ^ 

Since  the  reduced  problem  is  generally  nonlinear  and  since  COLSYS  solves 
nonlinear  problems  using  a  damped  Newton  method  we  have  to  supply  formulas  for 
evaluating  the  Jacobians  of  f,  Y,  0,  and  Y  with  respect  to  X.  We  do  this  by 
providing  analytical  formulas  for  these  Jacobians  that  neglect  the  influence 
of  the  derivatives  of  E,  L,  R,  and  G2.  This  procedure  has  not  failed  on  any 
if  our  examples;  however,  an  alternate  possibility  would  be  to  approximate  the 
Jacobians  by  finite  differences. 

C|K.  K.  O'Malley,  Jr.,  and  J.  E.  Flaherty,  "Numerical  Methods  For  Stiff  Systems 
of  Two-Point  Boundary  Value  Problems,"  to  appear. 

^B .  Childs,  M.  Scott,  J.  W.  Daniel,  E.  Denman,  and  P.  Nelson  (Eds.),  Codes  For 
Boundary-Value  Problems  in  Ordinary  Differential  Equations,  May  14-17,  197b, 
Lecture  Notes  in  Computer  Science,  No.  76,  Springer-Verlag ,  Berlin,  1979. 

7ij.  Ascher,  I.  Christiansen,  and  R.  0.  Russell,  "Collocation  Software  For 
Boundary  Value  OUE's,"  ACM  Trans.  Math.  Software,  7  (1981),  pp.  209-222. 


We  start  the  Newton  iteration  with  a  uniform  mash  and  the  default  initial 
guess  x(0)(t)  for  X(t)  that  is  provided  by  COLSYS  and  calculate  successive 
approximatlqms  x(p)(t)  until  convergence  is  attained.  At  each  iteration  step 
we  calculate  an  approximation  E(p)(t)  to  E(t)  for  t  ■  0  and  1  as  the  Schur 
decomposition  of  G2(x(P)(t) ,t ,0) .  In  the  examples  of  Section  4  we  used 
analytical  formulas  for  E  rather  than  the  numerical  procedures  of  Golub  and 
Wilkinson 3  or  Ruhe.^  Finally,  L^P)  and  R.(p)  are  obtained  using  Gaussian 
elimination  to  row  reduce  A2(X^ P)(0) ,0)E_(p)(0)  and  B2(x(p)( 1) ,0)E+(p)( 1) , 
respectively. 

When  the  above  procedure  converges  we  calculate  boundary  layer 
corrections  u(t)  and  \>(o),  for  a  given  value  of  c,  using  Eqs.  (12),  (16a), 
(19),  and  (22a),  and  add  these  to  the  reduced  solution  in  order  to  get  the 
0(e)  asymptotic  approximation  (Eq.  (3)).  For  moderately  small  values  of  e 
this  approximation  may  not  provide  a  sufficiently  accurate  representation  of 
the  solution  and,  in  this  case,  we  use  it  as  an  initial  guess  to  COLSYS  and 
solve  the  complete  problem  (Eq.  (1)).  Unfortunately,  this  procedure  will  fail 
unless  we  also  provide  COLSYS  with  an  initial  nonuniform  partition 

n  :■  (0  *  t0  <  ti  <  ...  <  tfl  *  1}  (23) 

that  is  appropriately  graded  within  the  boundary  layers.  We  seek  to  find  n  so 
that  the  pointwise  error  satisfies 

||e(ti)||  <  5(1  +  |  ju(ti)H)  ,  i  -  1,2 . N-l  (24) 

^G.  H.  Golub  and  J.  H.  Wilkinson,  "Ill-Conditioned  Eigensystems  and  the 
Computation  of  the  Jordon  Canonical  Form,”  SIAM  Review  18  (1976),  pp. 

578-619. 

^A.  Ruhe,  "An  Algorithm  for  Numerical  Determination  of  the  Structure  of  a 
General  Matrix,"  BIT,  10  (1970),  pp.  196-216. 
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where  5  ts  a  prescribed  tolerance,  u*"  :=  [x^,yT]  ,  e  is  the  difference  between 
u  and  its  collocation  approximation,  and 


||u(ti)||  :=  max  |uj(ti)|  (25) 

l<j<m+n 

We  have  based  our  condition  for  determining  it  on  a  potntwise  error  criteria 
since  this  seemed  to  work  better  in  practice  than  a  global  criteria.  This  is 
somewhat  surprising  since  COLSYS  uses  a  global  error  criteria  to  select  a 
mesh. 


We  assume  that  the  final  partition  selected  by  COLSYS  to  solve  the 
reduced  problem  satisfies  equation  (24)  outside  of  boundary  layer  regions  and 
we  seek  to  refine  it  within  the  boundary  layers.  We  further  assume  that 
derivatives  of  u  can  adequately  be  replaced  by  either  p(t)  or  v(o)  in  the  left 
or  right  boundary  layer,  respectively. 

This  problem  was  studied  by  Ascher  and  Weiss*  who  showed  that  Eq.  (24) 
could  he  approximately  satisfied  in  the  left  boundary  layer  by  choosing 
-’ubinterval  lengths  as 

c  5( 1+| ju(ti-i) | | )  */2k 


ti 


ti-i  =  (--)[- 


h(ti-i) 


(26) 


for  collocation  at  the  image  of  k  Gauss-Legend  re  points  per  subinterval .  Here 
o  is  a  numerical  constant  and  a_  is  the  magnitude  of  the  largest  diagonal 
•'lament  of  T_(X(0)).  A  similar  formula  can  be  obtained  for  selecting 
suh interval  lengths  In  the  right  boundary  layer. 


‘U.  Ascher  and  R.  Weiss,  "Collocation  For  Singular  Perturbation  Problems  1: 
First  Order  Systems  With  Constant  Coefficients,"  Technical  Report  81-2,  Dept. 
Comp.  Sci. ,  University  of  British  Columbia,  1981. 
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Starting  with  1  -  1  we  use  Eq.  (26)  to  generate  a  partition  until  we 
either  reach  t  ■  1/2  or  a  point  where  a  subinterval  length  selected  by  Eq. 

(26)  is  larger  than  that  used  by  COLSYS  to  solve  the  reduced  problem.  We  then 
repeat  the  procedure  in  the  right  boundary  layer. 

We  have  written  a  computer  code  called  SPCOL  that  implements  the 
algorithms  that  are  described  in  this  section;  thus,  it  (i)  uses  COLSYS  to 
solve  the  reduced  problem,  (ii)  calculates  and  adds  appropriate  boundary  layer 
corrections  to  the  reduced  problem,  and  (iii)  (optionally)  suggests  a  mesh 
that  can  be  used  by  COLSYS  to  solve  the  complete  problem. 

4 .  EXAMPLES 

In  order  to  appraise  the  performance  of  SPCOL  we  have  applied  it  to 

several  examples  involving  the  deformation  of  a  nonlinear  elastic  beam  which 

is  resting  on  a  nonlinear  elastic  foundation  and  is  subjected  to  the  combined 

action  of  a  horizontal  end  thrust  P  aud  a  lateral  load  p(x,t)  per  unit  length 

(cf.  Figure  1).  This  problem  is  discussed  and  analyzed  in  detail  in  Flaherty 

and  O'Malley^  and  herein  we  only  present  the  governing  equations,  which  in 

dimensionless  form  are 

•  •  • 

Xj  -  co 8  x  ,  X2  **  sin  X3  ,  X3  =  y\  (27a, b,c) 

eyi  »»  -y2  ,  ey2  “  (X^x2-p)  cos  x3  -  Tn  >  (27d,e) 

where 

T  =  sec  X3  +  ey2  tan  X3  (27i.i 


2j.  E.  Flaherty  and  R.  E.  O'Malley,  Jr.,  "Singularly-Perturbed  Boundary  Value 
Problems  for  Nonlinear  Systems,  Including  a  Challenging  Problem  For  a  Non¬ 
linear  Beam,"  Proceedings,  Conference  on  Singulare  Storungstheor ie  mit 
Anwendungen,  Oberwolfach,  1981. 
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The  slow  variables  (xj,x2)  and  X3  represent  the  Cartesian  coordinates  and  the 
tangent  angle  of  a  material  particle  on  the  centerline  of  the  beam  that  was  at 
the  Cartesian  location  (t,0)  in  the  undeformed  state.  The  fast  variables  yj 
and  y2  are  the  internal  bending  moment  and  transverse  shear  force, 
respectively  (cf.  Figure  1).  Finally,  the  small  parameter  is 

e2  «  El/ PL2  ,  (28) 

where  El  is  the  flexural  rigidity  and  L  is  the  length  of  the  beam;  thus,  our 
beam  is  much  stronger  in  extension  than  it  is  in  bending. 

This  example  does  not  precisely  fit  our  hypotheses  since  the  axial  force 
T  is  a  function  of  the  fast  variable  y2  and,  thus,  G2  also  depends  on  y. 
However,  our  theory  and  methods  will  still  apply  as  long  <>s  y  remains  bounded 
and  | x 3 1  <  ir/-2  as  £  +  0.  In  order  to  illustrate  the  diverse  behaviors  that 
can  occur  when  y  either  does  or  does  not  remain  bounded  as  e  +  0  we  present 
solutions  for  two  problems  both  having  A  =  p  =  1  and  which  differ  only  in 
their  boundary  conditions.  Some  additional  examples  are  presented  in  Flaherty 

and  O'Malley. 

In  our  first  example  we  take  the  boundary  conditions  as 

xi(0)  *  0  ,  -10x2(0)  +  y2(0)  *  0  ,  -X3(0)  +  10yi(0)  =  0 

(29) 

10x2(1)  +  y2( 1)  “  0  ,  10x3(1)  +  y 1 ( 1 )  ■  0 


2J.  E.  Flaherty  and  R.  E.  O'Malley,  Jr.,  "Singularly-Perturbed  Boundary  Value 
Problems  for  Nonlinear  Systems,  Including  a  Challenging  Problem  For  a  Non¬ 
linear  Beam,"  Proceedings,  Conference  on  Singulare  Storungstheorie  mlt 
Anwendungen ,  Oberwolfach,  1981. 

^R.  E.  O'Malley,  Jr.,  and  J.  E.  Flaherty,  "Numerical  Methods  For  Stiff  Systems 
of  Two-Point  Boundary  Value  Problems,"  to  appear. 
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These  supports  correspond  the  a  beam  that  is  almost  simply  supported  at  t  = 
and  almost  clamped  at  t  ■  1.  However,  perhaps  due  to  friction,  there  is  som 
coupling  between  lateral  and  rotational  effects  at  the  supports. 

As  we  shall  see,  y  remains  bounded  in  this  example  so  our  methods  are 
applicable.  The  orthogonal  matrix 


E(x(0)) 


(l+a2)-l/2 


1 


-a 


(  It);.  ) 


a  1 


where 


a2  ■  sec  X3(0) 


(  lot-,) 


reduces 


to  the  Schur  form  given  by  equation  (7)  at  t  -  0  and  ET  will  reduce 
G2(x(1),1,0)  to  the  form  given  by  Eq.  (17)  at  t  *  1. 

We  solved  this  problem  in  two  ways:  (i)  using  COLSYS  to  solve  the 
complete  problem  (Eqs.  (27)  and  (29))  with  continuation  from  a  large  to  a 
small  value  of  e  and  (ii)  using  our  code  SPCOL  to  compute  an  initial 
asymptotic  approximation  and  to  recommend  a  nonuniform  mesh  and  using  this 
with  COLSYS  to  calculate  an  improved  solution.  All  calculations  were 
performed  in  double  precision  on  an  IBM  3033  computer,  used  two  collocation 
points  per  subinterval,  and  set  the  error  tolerance  5  (cf.  Eq.  (24))  at  l>'-b 
for  slow  variables  and  10“3  for  fast  variables. 

Our  results  for  the  normalized  CP  times  and  the  number  of  subintervals 
(NSUB)  that  are  either  used  by  COLSYS  or  recommended  by  SPCOL  are  shown  in 


G2(x(0) ,0,0) 


-u2 
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Tables  I  and  II  for  continuation  in  e  and  our  methods,  respectively.  Differ¬ 
ences  between  our  initial  asymptotic  approximation  and  the  final  solution 
obtained  by  COLSYS  are  shown  for  X3  and  y2  at  t  *  0  and  1  in  Table  III.  We 
see  that  the  differences  decrease  like  0(e)  as  expected.  Differences  that  are 
recorded  as  zero  are  less  than  10“8.  Finally,  we  exhibit  solutions  for  X2, 
x3»  yi»  at»d  Y2  *n  Figure  2. 

TABLE  I.  NONLINEAR  ELASTICALLY  SUPPORTED  BEAM.  NUMBER  OF  SUBINTERVALS  (NSUB) 
AND  CP  TIMES  USED  TO  SOLVE  THE  PROBLEM  BY  COLSYS  WITH  CONTINUATION 
IN  e.  THE  TOTAL  CP  IS  THE  ACCUMULATED  TIME  FOR  THE  e  SEQUENCE. 


NSUB 

CP 

Total  CP 

80 

8.0 

8.0 

78 

9.0 

17.0 

78 

19.5 

36.5 

156 

44.5 

81.0 

100 

19.0 

100.0 

TABLE  II.  NONLINEAR  ELASTICALLY  SUPPORTED  BEAM.  NUMBER  OF  SUBINTERVALS 

(NSUB)  AND  CP  TIMES  TO  SOLVE  THE  PROBLEM  BY  SPCOL  AND  OBTAIN  AN 
IMPROVEMENT  BY  COLSYS.  THE  CP  TIMES  FOR  SPCOL  INCLUDE  THE  TIME  TO 
CALCULATE  THE  REDUCED  SOLUTION  WHICH  WAS  4.8  TIME  UNITS. 

CORRECTION  NO.  1  USES  THE  MESH  THAT  WAS  RECOMMENDED  BY  SPCOL. 
CORRECTION  NO.  2  USES  A  MESH  THAT  IS  TWICE  AS  COARSE.  THE  TOTAL  CP 
IS  THE  SUM  OF  THE  TIMES  FOR  THE  SPCOL  AND  COLSYS  SOLUTIONS. 


SPCOL 


Rec.  No. 
of  NSUB 


COLSYS 

Correction  No.  1 


NSUB  CP  Total  CP 


10-1 

40 

10-2 

45 

10-4 

54 

10~8 

55 

COLSYS 

Correction  No.  2 

NSUB 

CP 

Total  CP 

80 

12.1 

16.9 

78 

8.1 

12.9 

66 

9.2 

14.1 

Failed 

TABLE  III.  NONLINEAR  ELASTICALLY  SUPPORTED  BEAM.  DIFFERENCES  BETWEEN  SPCOL 
AND  COLSYS  SOLUTION &,  WHERE  A(  )  |(  )SPCOl  -  (  >COLSYsl 


c 

Ax3(0) 

Ay2(0) 

Ax3( 1) 

Ay2( 1) 

io-1 

3.3xl0_1 

5. lxlO-2 

6.8xl0-1 

3.6xl0_1 

10"2 

2.8xlQ-2 

6.6xlO“3 

6.  lxlO-2 

3.9xl0-2 

10”4 

2.7xl0“4 

6.8xlO"5 

6.  lxlO-4 

3.9xl0*4 

10~8 

0 

1. 3xl0“7 

0 

0 

The  results  reported  In  Tables  I  and  II  need  some  additional  explanation. 
The  number  of  subintervals  and  CP  times  used  with  continuation  depended  heav¬ 
ily  on  the  e  sequence  that  was  used.  The  results  in  Table  I  are  about  the 
best  insofar  as  they  gave  the  smallest  total  CP  time  for  the  sequence.  In 
addition,  COLSYS  relies  on  the  difference  between  solutions  that  are  computed 
on  two  different  partitions  in  order  to  estimate  local  errors.  Thus,  at  a 
minimum,  COLSYS  would  always  double  our  suggested  mesh.  This  is  apparent  in 
the  results  listed  under  the  heading  of  “COLSYS  Correction  No.  1"  in  Table  II. 
In  some  sense  these  results  are  encouraging  insofar  as  they  indicate  that  our 
mesh  selection  strategy  is  doing  about  as  well  as  it  can,  at  least  for 
e  <  10-2  .  However,  it  seems  that  fewer  points  should  be  necessary,  so  we 
tried  giving  COLSYS  an  initial  mesh  that  consisted  of  every  other  point  of  our 
recommended  mesh.  This  is  clearly  a  risky  strategy  since  collocation  at  the 
Gauss-Legendre  points  is  known  to  be  unstable  unless  the  mesh  is  sufficiently 
fine  in  the  boundary  layers  (cf.  Ascher  and  Weiss^).  Our  results  using  this 

*U.  Ascher  and  R.  Weiss,  "Collocation  For  Singular  Perturbation  Problems  I: 

First  Order  Systems  With  Constant  Coefficients,"  Technical  Report  81-2,  Dept. 

Comp.  Scl.,  University  of  British  Columbia,  1981. 


are  reported  under  the  heading  of  "COLSYS  Correction  No.  2"  in  Table  II.  Some 
improvement  is  noted  for  e  >  10— **  ;  however,  COLSYS  failed  to  find  a  solution 


(within  our  prescribed  limitations)  when  e  -  10-8. 

In  our  second  example  we  use  the  boundary  conditions 

xi(0)  =>  0  ,  .-X2(°)  +  ey2<0)  **  0  ,  -X3(0)  +  e2yi(0)  *  0 


(32) 

X2(l)  +  ey2(0  ■  0  ,  X3O)  +  e2yi(l)  =  0 

If  e  were  set  to  zero  then  these  boundary  conditions  would  correspond  to 
clamped  supports  at  t  «  0  and  1.  Since  the  limiting  boundary  conditions  only 
involve  the  slow  variables  and  since  the  slow  vector  x  cannot  generally 
satisfy  all  of  them  as  e  +  0  we  would  expect  the  solution  to  have  boundary 
layers  in  these  components.  This  in  turn  will  force  the  fast  vector  y  to 
become  unbounded  like  0(l/e)  at  the  endpoints!  Thus,  this  problem  does  not 
have  an  asymptotic  expansion  having  the  form  of  Eq.  (3);  however,  an 
appropriate  asymptotic  representation  of  a  solution  has  been  obtained  by 
Flaherty  and  O'Malley.2  We  shall  not  repeat  those  results  here,  but  in  order 
to  emphasize  the  diverse  behavior  that  can  occur  in  nonlinear  singularly- 
p^rfnrbed  problems,  we  present  solutions  for  X2,  X3,  eyj ,  and  ey2  in  Figure  3 
These  solutions  were  computed  using  COLSYS  with  continuation  in  e. 


5.  DISCUSSION 

We  have  obtained  asymptotic  approximations  for  a  restricted  class  of 
nonlinear  singularly-perturbed  boundary  value  problems  and  have  shown  how  to 


2J.  E.  Flaherty  and  R.  E.  O'Malley,  Jr.,  "Singularly-Perturbed  Boundary  Value 
Problems  For  Nonlinear  Systems,  Including  a  Challenging  Problem  For  a  Non¬ 
linear  Beam,"  Proceedings,  Conference  on  Singulare  Storungstheorle  mit 
Anwendungen,  Oberwolfach,  1981. 


construct  them  numerically  and  use  them  to  suggest  a  nonuniform  mesh  that  may 
be  used  as  input  to  a  two-point  boundary  value  code  in  order  to  calculate 
Improved  solutions.  Clearly  this  approach  offers  some  advantages  over  the 
more  standard  technique  of  continuation  in  e  steps;  however,  the  picture  is 
far  from  clear  and  several  questions  still  remain  as  to  how  best  to  use 
asymptotic  analysis  in  conjunction  with  numerical  analysis. 

As  we  have  shown  in  our  second  example  of  Section  4,  very  diverse 
behavior  in  the  solution  of  singularly-perturbed  problems  can  result  from 
seemingly  minor  changes  in  boundary  conditions.  Some  phenomena  cannot  easily 
be  predicted,  so  perhaps  a  sensible  course  to  follow  is  to  use  asymptotic  and 
numerical  methods  in  tandem.  For  example,  a  rough  numerical  solution  could  be 

v 

obtained  for  several  values  of  e  which  could  then  be  used  to  suggest  the  form 
of  an  asymptotic  solution.  The  asymptotic  approximation  could  then  be  used  to 
refine  the  numerical  solution,  and  so  on.  It  is  also  possible  that  singular 
perturbation  theory  could  be  used  to  construct  special  methods  that  are 
appropriate  for  specific  problems  as  e.g.,  in  Flaherty  and  Mathon®  and 
Ascher  and  Weiss. ^ 

Throughout  our  discussion  we  have  ignored  the  question  of  uniqueness.  In 
general,  multiple  solutions  can  be  expected  and  they  must  be  coped  with 

*U.  Ascher  and  R.  Weiss,  "Collocation  For  Singular  Perturbation  Problems  Is 
First  Order  Systems  With  Constant  Coefficients,"  Technical  Report  81-2,  Dept. 
Comp.  Scl.,  University  of  British  Columbia,  1981. 

8j.  E.  Flaherty  and  W.  Mathon,  "Collocation  With  Polynomial  and  Tension 
Splines  For  Singularly-Perturbed  Boundary  Value  Problems,"  SIAM  J.  Scl.  Stat. 
Comput. ,  J^,  (1980),  pp.  260-289. 


numerically.  In  Reference  9  we  showed  how  asymptotic  methods  may  be  used  to 
distinguish  the  different  solutions  and  to  provide  initial  guesses  for  a 
two-point  boundary  value  code. 


9J.  E.  Flaherty  and  R.  E.  O’Malley,  Jr.,  “On  the  Numerical  Integration  of  Two- 
Point  Value  Problems  For  Stiff  Systems  of  Ordinary  Differential  Equations," 
Boundary  and  Interior  Layers  -  Computational  and  Asymptotic  Methods,  J.  J.  H. 
Miller’,  Editor,  Boole  Press,  Dublin,  1980,  pp.  93-102. 
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DIRECTOR,  OPERATIONS  DIRECTORATE  1 

DIRECTOR,  PROCUREMENT  DIRECTORATE  1 

DIRECTOR,  PRODUCT  ASSURANCE  DIRECTORATE  1 


NOTE:  PLEASE  NOTIFY  DIRECTOR,  BENET  WEAPONS  LABORATORY,  ATTN:  DRDAR-LCB-TL, 
OF  ANY  REQUIRED  CHANGES. 


TECHNICAL  REPORT  EXTERNAL  DISTRIBUTION  LIST 


NO.  OF  NO.  OF 

COPIES  COPIES 


ASST  SEC  OF  THE  ARMY 

RESEARCH  &  DEVELOPMENT 

ATTN:  DEP  FOR  SCI  &  TEQH  1 

THE  PENTAGON 

WASHINGTON,  D.C.  203L5 

COMMANDER 

DEFENSE  TECHNICAL  INFO  CENTER 
ATTN:  DTIC-DDA  12 

CAMERON  STATION 
ALEXANDRIA,  VA  22314 

COMMANDER 

US  ARMY  MAT  DEV  &  READ  COMD 
ATTN :  DRCDE-SC 
5001  EISENHOWER  AVE 
ALEXANDRIA,  VA  22333 

COMMANDER 
US  ARMY  ARRADCOM 
ATTN:  DRDAR-LC 

DRDAR-LCA  (PLASTICS  TECH 
EVAL  CEN) 

DRDAR-LCE 

DRDAR-LCM  (BLDG  321) 

DRDAR-LCS 
DRDAR-LCU 
DRDAR-LCW 
DRDAR-TSS  (STINFO) 

DOVER,  NJ  07801 

DIRECTOR 

US  ARMY  BALLISTIC  RESEARCH  LABORATORY 
ATTN:  DRDAR-TSB-S  (STINFO) 

ABERDEEN  PROVING  GROUND,  MD  21005 

COMMANDER 
US  ARMY  ARRCOM 
ATTN:  DRSAR-LEP-L 

ROCK  ISLAND  ARSENAL 
ROCK  ISLAND,  IL  61299 


COMMANDER 

ROCK  ISLAND  ARSENAL 

ATTN:  SARRI-ENM  (MAT  SCI  DIV)  1 

ROCK  ISLAND,  IL  61299 

DIRECTOR 

US  ARMY  INDUSTRIAL  BASE  ENG  ACT 
ATTN:  DRXIB-M  1 

ROCK  ISLAND,  IL  61299 

COMMANDER 

US  ARMY  TANK-AUTMV  R&D  COMD 
ATTN:  TECH  LIB  -  DRSTA-TSL  1 

WARREN,  MI  48090 

COMMANDER 

US  ARMY  TANK-AUTMV  COMD 
ATTN:  DRSTA-RC  1 

WARREN,  MI  48090 

1  COMMANDER 

1  US  MILITARY  ACADEMY 

ATTN:  CHM,  MECH  ENGR  DEPT  1 

1  WEST  POINT,  NY  10996 

1 

1  US  ARMY  MISSILE  COMD 

1  REDSTONE  SCIENTIFIC  INFO  CEN 

1  ATTN:  DOCUMENTS  SECT,  BLDG  4484  2 

2  REDSTONE  ARSENAL,  AL  35898 

COMMANDER 

US  ARMY  FGN  SCIENCE  &  TECH  CEN 
ATTN:  DRXST-SD  1 

1  220  7TH  STREET,  N.E. 

CHARLOTTESVILLE,  VA  22901 

COMMANDER 

US  ARMY  MATERIALS  &  MECHANICS 
1  RESEARCH  CENTER  2 

ATTN:  TECH  LIB  -  DRXMR-PL 

WATERTOWN,  MA  02172 


NOTE:  PLEASE  NOTIFY  COMMANDER,  ARRADCOM,  ATTN:  BENET  WEAPONS  LABORATORY, 

DRDAR-LCB-TL,  WATERVLIET  ARSENAL,  WATERVLIET,  NY  12189,  OF  ANY 
REQUIRED  CHANGES. 


TECHNICAL  REPORT  EXTERNAL  DISTRIBUTION  LIST  (CONT'D) 


NO.  OF 
COPIES 

COMMANDER 

US  ARMY  RESEARCH  OFFICE 

ATTN:  CHIEF,  IPO  1 

P.O.  BOX  12211 

RESEARCH  TRIANGLE  PARK,  NC  27709 
COMMANDER 

US  ARMY  HARRY  DIAMOND  LAB 

ATTN:  TECH  LIB  1 

2800  POWDER  MILL  ROAD 

ADELPHIA,  MD  20783 

COMMANDER 

NAVAL  SURFACE  WEAPONS  CEN 
ATTN:  TECHNICAL  LIBRARY  .  1 

CODE  X212 

DAHLGREN,  VA  22448 


NO.  OF 
COPIES 

DIRECTOR 

US  NAVAL  RESEARCH  LAB 
ATTN:  DIR,  MECH  DIV  1 

CODE  26-27  (DOC  LIB)  1 

WASHINGTON,  D.C.  20375 

METALS  &  CERAMICS  INFO  CEN 
BATTELLE  COLUMBUS  LAB 
505  KING  AVE  1 

COLUMBUS,  OH  43201 

MATERIEL  SYSTEMS  ANALYSIS  ACTV 
ATTN:  DRXSY-MP 

ABERDEEN  PROVING  GROUND  L 

MARYLAND  21005 


NOTE:  PLEASE  NOTIFY  COMMANDER,  ARRADCOM,  ATTN:  BENET  WEAPONS  LABORATORY, 
DRDAR-LCB-TL,  WATERVLIET  ARSENAL,  WATERVLIET,  NY  12189,  OF  ANY 
REQUIRED  CHANGES. 


